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ABSTRACT 


Numerical  methods  for  solving  the  heat  equation  via  potential  theory  have  been  ham¬ 
pered  by  the  high  cost  of  evaluating  heat  potentials.  When  M  points  are  used  in  the 
discretization  of  the  boundary  and  N  time  steps  are  computed,  an  amount  of  work  of  the 
order  0(N2M 2)  has  traditionally  been  required.  In  this  paper,  we  present  an  algorithm 
which  requires  an  amount  of  work  of  the  order  O(NM),  and  we  observe  speedups  of  five 
orders  of  magnitude  for  large-scale  problems.  Thus,  the  method  makes  it  possible  to  solve 


1  Introduction 

A  classical  approach  to  the  solution  of  the  heat  equation 


Ut  =  A  U 


in  a  space-time  domain  Cl?  =  Y\t=o  (see  Fig.  1)  is  through  the  use  of  heat  potentials 
[3,7].  Given  zero  initial  conditions,  one  seeks  a  representation  of  U  as  a  single  layer  heat 
potential 

Sfi(x,t)  =  f  f  K(x,x' ,t  —  f'Wx' ,  f')  dx'  dt'  (1) 

Jo  Jr(t') 

or  a  double  layer  heat  potential 

Dn(x,t)=[  f  -t')ii{x! ,t')dx!  dt'  ,  (2) 

Jo  Vr(t')  on' 

where  if  is  a  fundamental  solution  of  the  heat  equation  in  some  region  containing  f2,  n' 
denotes  the  unit  outward  normal  to  T(t')  =  df l(t')  at  x',  and  fj.  is  a  surface  density  defined 
on  Tt  =  nf=o 


Figure  l 

A  domain  Qj  contained  in  space-time  Rn  X  R.  f 2(f)  is  the  cross-section  of  Qt  at  time  t-  Tj  is  the 
lateral  boundary  of  0.T-  Its  cross-section  at  time  t  is  given  by  T(£)  = 
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For  example,  the  initial-Neumann  problem  where  • 

U(x,  0)  =  0  in  S7(0) 

dU  r 

a^  =  9  on  Tt 

is  reduced  to  the  Volterra  integral  equation  of  the  second  kin  i 

\v(x,t)  +  Dfi(x,t)  =  g{x,t)  on  Ft 

by  means  of  the  representation  U  =  Sp.  Such  an  approach  is  used  for  numerical  solution 
of  the  heat  equation  in  [4,6]. 

Other  applications  of  heat  potentials  include  the  study  of  crystal  growth  and  unstable 
solidification  which  can  be  modelled  by  the  integral  equation 

eC  +  V  +  U  +  SV  =  0, 

where  the  unknown  V  is  the  normal  velocity  of  the  solid-liquid  boundary  T(t),  C  is  the 
curvature  of  T(t),  e  is  a  positive  constant,  and  U  is  a  given  temperature  field  incorporating  f 

the  initial  and  boundary  conditions.  Recent  numerical  methods  for  unstable  solidification 
have  used  this  formulation  [5,8], 

However,  numerical  methods  based  on  heat  potentials  have  been  crippled  by  their 
history-dependence:  both  solving  the  integral  equation  and  evaluating  the  potential  rep¬ 
resentation  require  more  and  more  work  as  time  proceeds.  Consider,  for  example,  the 
task  of  calculating  S/i  at  a  sequence  of  time  levels  t  =  At,  2At, . . . ,  NAt.  At  the  nth 
level,  we  must  sum  over  n  previous  levels;  therefore,  the  total  work  is  0(N2).  As  for  the 
spatial  variables,  if  we  are  given  M  points  Xj,  j  =  1,2, ... ,M  in  the  discretization  of  the 
boundary,  the  evaluation  of  Sfi  requires  0(M2)  work  per  time  step.  Thus,  solving  the  heat 
equation  up  to  a  fixed  time  T  by  an  integral  equation  method  would  seem  to  require  at 
least  0(N2M2)  work.  This  high  cost  has  prevented  integral  equation  methods  from  being 
used  in  practice. 

In  this  paper,  we  develop  an  algorithm  for  the  rapid  evaluation  of  heat  potentials.  This 
algorithm  requires  only  O(MN)  work  to  evaluate  or  Du  at  M  points  on  the  boundary 
at  each  of  N  time  levels.  Since  there  are  M N  data  points  and  MN  values  to  be  computed, 
this  is  asymptotically  optimal.  The  basic  idea  of  the  algorithm  is  that  the  potential  can  be 
split  into  two  components,  one  representing  the  effects  of  the  source  n  over  distant  time 
(the  history  part)  and  one  representing  the  effects  of  /i  over  recent  time  (the  local  part). 

The  history  part  is  smooth  and  can  be  well  approximated  by  only  a  few  Fourier  modes. 

The  local  part,  on  the  other  hand,  can  be  well  approximated  by  Taylor  expansion,  using 

the  singularity  structure  of  the  fundamental  solution.  * 

The  outline  of  the  paper  is  as  follows;  Section  2  describes  the  fundamental  solution  for 
a  box,  Section  3  describes  the  fast  algorithm  itself,  and  Section  4  presents  some  numerical 
results.  We  state  our  conclusions  in  Section  5. 
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The  algorithm  is  currently  being  applied  to  crystal  growth  and  unstable  solidification 
by  Sethian  and  Strain  [9].  In  a  subsequent  paper,  we  will  describe  a  different  local  approx¬ 
imation  and  use  the  fast  algorithm  for  the  numerical  solution  of  the  heat  equation. 

A  detailed  description  of  the  algorithm  is  presented  only  for  the  single  layer  potential 
Sp\  the  modifications  necessary  for  evaluating  Dp  or  volume  potentials  are  straightforward. 
Furthermore,  we  consider  only  potentials  formed  with  the  “box”  kernel  defined  below;  one 
advantage  of  the  integral  equation  approach  is  that  the  solution  in  f It  does  not  depend  on 
which  kernel  is  used. 


2  The  Fundamental  Solution  in  a  Box 

In  this  section,  we  obtain  complementary  representations  for  the  fundamental  solution  K 
of  the  heat  equation  in  a  box  B  =  [0,  l]n  with  homogeneous  Dirichlet  boundary  conditions. 
A  Fourier  series  calculation  [2]  shows  that 

A'(x,x',t)  =  2n  ^  e-!f2lkl2*  sm(Tr k{X i) sin(Tc kix'i)  ,  (3) 

«  keNn  «=i 

where  N  denotes  the  positive  integers,  k  —  (fci, kn )  and  x  =  (aq, ...,  xn). 

On  the  other  hand,  the  method  of  images  (see  [10])  can  be  used  to  show  that 

A'(x,x',()  =  (4irf)'"'/2  £  E  ■■■  E  ,  (4) 

keZ"  i7i=±i  <7„=±i 

where  a  •  x'  =  (<7i  x\, . . . ,  an  x'n).  This  expression  can  also  be  derived  by  the  Poisson 
summation  formula  [2]. 

The  equality  of  (3)  and  (4)  is  one  of  the  foundations  of  our  algorithm;  both  sums 
converge  exponentially  fast,  but  in  different  regions  of  time.  Indeed,  the  error  in  using  pn 
terms  of  (3)  is  of  the  order  0(e~p2ir2t )  as  p2t  — >  oo,  whereas  the  error  in  using  2n(2 p  +  l)n 
terms  of  (4)  is  of  the  order  0(e-p2/{)  as  p2 Jt  — >  oo.  Thus,  we  truncate  the  Fourier  series 
(3)  to  evaluate  K  for  large  t  and  the  sum  of  Gaussians  (4)  to  evaluate  K  for  small  t. 


3  The  Fast  Algorithm 

The  fast  algorithm  will  be  explained  in  the  simplest  context,  namely  that  of  evaluating  the 

#  single  layer  heat  potential 

Sp(x,t)=  f  f  K(x,x’ ,t  -  t')p.{x' ,t')dx'  dt'  (5) 

*  Jo  Jr(t') 

in  two  space  dimensions.  Here  T(t)  is  a  family  of  boundary  curves  lying  in  the  unit  box 
B  =  [0,  l]2  and  K  is  the  fundamental  solution  of  the  heat  equation  in  B  with  homogeneous 
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Dirichlet  boundary  conditions: 


A(x,x',t)  =  4  ^2  e  ,r2'k'2t  sin(7rfciX1)sin(7rA:1a:,1)sin(7rA:2a:2)sin(7rfc2a:/2)  (6) 

keN* 


=  FT  E  E  *1*2  e-»*— 2k"2/4t  - 

keZ2  <u=±i 

(7) 

We  begin  by  splitting  the  time  integral  in  (5)  at  time  t  —  8,  with  8 
be  determined  later.  Thus,  we  write  S/i  =  Spfi  +  Spfi,  where 

a  small  parameter 

5^(x,t)  =  f  f  A'(x,  x' ,t  —  t'Wx' ,t')  dx'  dt1  , 
Jt-s  Jr(t') 

(8) 

SpiJ.{x,t)  =  I  [  A’(x,x' ,t  —  t'Wx' ,t')  dx'  dt1 

Jo  Jr(t') 

(9) 

The  subscripts  L  and  F  refer  to  the  local  and  Fourier  parts,  respectively. 


3.1  Fast  evaluation  of  Spjt 

First  consider  the  component  Spfi,  which  contains  the  history-dependence  of  the  potential.  * 

After  replacing  K  with  its  Fourier  expansion  (6),  Spu  becomes  a  Fourier  series 

oo  oo 

SFn(x,t)  =  Ck(t,  6)  sin(7rfc1x1)sin(7rfc2X2)  ,  (10) 

fcl=l fc2=l 

with  coefficients 

Ck(*»£)  =  4  f  e— ^-3|k|2(f— t')  f  sin(7rA:1x,1)sin(7r^2^2)/x(x,T,)<^x,  dt' .  (11) 

Jo  Jr(t') 

This  representation,  by  itself,  does  not  eliminate  the  problem  of  history-dependence,  be¬ 
cause  each  of  the  Fourier  coefficients  Ck  at  time  t  is  obtained  by  integrating  over  all 
previous  history.  However,  Ck(t,  8)  can  be  computed  from  Ck(t  —  At,  6)  recursively.  That 
is,  by  separating  the  final  time  interval  At  from  the  rest,  we  get 

Ck(t,8)  =  e-*'WAtCk(t-At,8)  + 

4  /  e-’r2lkl2(<~t )  /  sin(n  kix\)  sin(  n  k2x'2)  t')  dx' dt'.  (12) 

Jt-At-s  Jr(v) 

Each  coefficient  Ck  can  be  updated  with  constant  work  per  time  step,  rather  than  recom-  4 

puted  from  t  —  0;  history-dependence  is  effectively  eliminated. 

Remark  3.1  This  elimination  of  history  actually  applies  much  more  generally.  Let  an  operator  A  I 

generate  a  semigroup  e‘A,  and  consider  the  problem  of  evaluating  the  Duhamel  integral 

F(t)  =  f*  e(£-t,)i4/(t')  dt' 

Jo 
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at  a  sequence  of  time  levels  t  =  At,  2At, . . . ,  N  At.  Redoing  the  time  integral  at  each  step  costs  0(N2)  work, 
but  we  can  compute  F(t  +  At)  recursively  from  F(t ): 


F(t  +  A<) 


/•f+At 

Jo 


> Af(t’)dt ' 


(13) 

(14) 


This  costs  only  O(N)  work  up  to  time  N At.  In  the  present  paper,  A  is  the  Laplacian  A  on  the  box  B  with 
Dirichlet  boundary  conditions  on  dB,  and  f(t')  is  a  measure  concentrated  on  T(t')  with  density  p. 


Another  feature  of  the  Fourier  series  representation  (10)  for  Spp  is  that  it  allows  us  to 
take  advantage  of  the  smoothing  effect  of  the  heat  operator.  Higher  modes  are  damped 
exponentially,  so  that  the  Fourier  series  representation  (10)  of  the  kernel  can  be  truncated 
after  p2  terms  with  an  exponentially  small  error. 

Lemma  3.1  Let  t  —  t‘  >  8  and  let  Ep  be  the  error  in  truncating  the  series  expansion  (6) 
after  p2  terms 

Ep  =  lA'Xx,  x',t  —  O  —  4  E  E  sin(7rfc,ii)  sin(7rfc,-x(-)|  . 

k\  =1  At2  =  1 


Then 


Proof: 


i=i 


Ep< 


e-2**(p+l)*« 
7T  6 


(15) 


oo  oo 


Ep  =  4|  E  E  e  7r2'k'2^  1 ^  JJ  sin(7rA;,x,)sin(7rAr,x()| 

fcl=P+lfc2=P+l  1  =  1 

OO  OC  /  CO  ^  ^ 

4  E  E  =  ( 2  E  ^kls 

*1=P+1  fc2=pH-l 


< 


k\  =p+I 


<  2e~*2(p+1^s  E  e~’r2]2s 


< 


3=0 

^2e-ir2(p+1)2fi  y°°  e-,r2|5'x2  \ 

( 2e~*2{p+l)2s  •  1  ^  ^  =  _Le-2>r2(p+l)J5 

\  2  \Ztt2S  ) 


(16) 


(17) 


An  error  bound  for  truncation  of  the  Fourier  series  representation  of  Spp  now  follows 
immediately. 


5 


Lemma  3.2  Let  Ep(p )  be  the  error  in  truncating  the  series  expansion  (10)  after  p2  terms 

p  p 

Ef(p)  =  |5f//(x, <)  —  ]T  ^  Ck(M)sin(7rfciXi)sin(7rfc2x2)|  • 

fci=i fc2=i 

Then 

EF(p)  <  Ellkke-a^tp-n)^  ,  (13) 

7 TO 

where  |Tt|  m  t/ie  area  o/Tt  and  |/i|oo  m  the  maximum  of  \p\  over  Ft- 
We  now  define  the  updates  Uk(t,  At,  6)  by 

Uk(t,  At,  8)  —  4  f  e— jt  |lc|2(t— t )  f  sin(7rA;1Xl)sin(7rA:2x2)^(x,,f,)  dx'  dt'  .  (19) 

Jt-s-At  xr(t') 

As  the  calculation  of  Sp  proceeds  in  time,  the  recursion  relation  (12)  provides  us  with 
a  means  for  updating  the  representation  of  Spp  at  a  total  cost  which  grows  only  linearly 
with  the  number  of  time  steps,  rather  than  quadratically: 

Ck(nAt,  S)  =  e-*W4‘ck((?i  -  1)A*,  S)  +  Uk(nAt,  At,  6)  .  (20) 

We  therefore  avoid  both  the  computational  cost  and  excessive  storage  required  by  the 
direct  evaluation  of  the  integral  Spp. 

We  still  need  to  construct  space  and  time  quadratures  for  evaluating  the  updates  Uk 
in  equation  (19).  First,  consider  the  calculation  of  the  trigonometric  moments  of  p 

yVdk(f,,^)=  /  sin^fcjx^sin^kjx^^xV^dx' , 

Jr(t') 

assuming  p  to  be  known  at  M  equidistant  points  on  r(f/).  The  trapezoidal  rule  for  smooth 
periodic  functions  converges  superalgebraically,  so  it  would  be  natural  to  use  it  to  evaluate 
Adk.  Using  all  M  points  to  integrate  each  of  the  first  p2  moments  at  each  time  step  leads 
to  a  nonoptimal  method,  however,  because  p  must  increase  as  8  — >  0  and  M  — +  oo,  to 
ensure  that  the  Fourier  series  truncation  error  vanishes  (see  equation  (18)  above).  But 
each  moment  Mk(t,6)  involves  integration  (over  a  smooth  curve)  of  sine  functions  with 
wave  numbers  kt  <  p.  Since  accurate  integration  of  such  an  oscillatory  function  requires 
only  a  fixed  number  of  quadrature  points  per  wavelength,  we  need  only  use  O(p)  of  the 
given  points  in  the  integration  scheme.  (The  constant  in  0(p )  will  depend,  of  course,  on 
the  smoothness  of  the  boundary  T  and  the  density  p.)  Hence,  all  of  the  p2  coefficients 
Ck(nAt,8)  can  be  updated  at  a  total  cost  of  0(p3)  work.  We  want  0(M )  work  per  time 
step,  so  we  choose  p  =  M 1^3.  In  the  error  estimate  (18),  we  then  have 

EF(p)  <  .  e-2*3Mh 

7 t8 
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By  choosing  8  =  M  3+e,  0  <  e  <  |,  vve  have 


Ef(p)  < 


Ml-*  -  |rT||^| 


CO  e~27 r2Me 


7T 


(21) 


which  is  decaying  superalgebraically  with  M.  For  example,  with  e  =  4/15  and  M  =  10, 
this  error  bound  is  already  less  than  10-15. 

Next  we  must  carry  out  the  time  integral  in  the  update.  Standard  approaches  like  the 
trapezoidal  rule  are  not  uniformly  second  order  accurate  as  N ,  M  — >  oo  and  <5  — ►  0,  because 
of  the  singularity  of  the  heat  kernel  at  t'  =  t,  which  lies  a  distance  8  away  from  the  endpoint 
of  the  interval  of  integration.  However,  we  can  evaluate  the  integral  of  an  exponential  times 
a  polynomial  exactly,  which  suggests  the  construction  of  a  product  integration  rule.  Thus, 
we  construct  weights  Wx  such  that  the  rule 


fl  &  e-*W dt'  =  W0g{t  —  At  —  6)  +  Wxg(t  -  8) 
Jt-At-S 


(22) 


is  exact  whenever  g  is  linear.  (Higher  order  rules  are  equally  easy  to  construct,  but  for 
simplicity,  we  are  seeking  a  globally  second  order  method.)  Some  simple  algebra  gives 


Wo  =  - 
Wi  =  e 


-1  +  2 


At  e_1r2|k|2s 


-1-2 


Ate-**WS 


(23) 

(24) 


where  2  =  7r2|k|2Al 

By  evaluating  the  updates  in  the  above  manner  (the  trapezoidal  rule  in  space  and 
product  integration  in  time),  it  is  clear  that  the  error  incurred  in  the  calculation  of  the 
Fourier  coefficients  Ck(t,<5)  is  0(At2)  plus  a  term  which  is  decaying  superalgebraically 
with  As  =  \T(t)\/M.  The  net  work  required  is  O(M)  per  time  step.  Now,  given  the  values 
Ck(t,5),  it  remains  only  to  evaluate  the  truncated  series 

p  p 

SFn(x,t )  =  iL,  Ck(<,^)sin(7rfcir1)sin(7rfc2a:2) 

fcl=l  &j=l 


at  the  M  points  Xj  given  on  T(t).  Direct  evaluation  of  the  series  would  require  0(Mp 2) 
work  which  would  preclude  optimality.  However,  SFp  contains  only  information  with 
wave  numbers  kt  <  p.  It  suffices,  therefore,  to  evaluate  it  directly  at  O(p)  of  the  given 
points  on  the  one-dimensional  set  T(t).  The  values  of  SFp  can  then  be  reconstructed 
at  the  rest  of  the  points  by  high-order  local  interpolation  along  the  curve.  Using  the 
preceding  strategy,  the  total  amount  of  work  required  is  of  the  form  0(p3)  +  O(M)  to 
achieve  some  fixed  interpolation  error,  say  second  order  in  As.  (However,  since  the  points 
on  the  curve  are  equispaced  in  arclength,  one  can  actually  make  the  interpolation  error 
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decay  superalgebraically,  by  means  of  an  FFT  in  arclength.  The  net  computational  cost 
would  then  be  of  the  form  0(p3)  +  0(M  log  M).) 

In  summary,  we  can  evaluate  Sf/.i  in  0(M)  work  per  time  step,  with  a  constant  inde¬ 
pendent  of  iV,  and  with  an  error  of  the  form 

0(At2  +  As2  +  M2^e-2^Me)  , 

when  8  =  0(M~2/3+( )  and  second  order  interpolation  is  used  in  arclength.  Since  the  final 
term  is  superalgebraic,  the  error  in  Spfi  is  second  order.  Higher  order  error  estimates  can 
be  obtained  by  using  higher  order  product  integration  and  higher  order  interpolation  in 
arclength. 


3.2  Fast  evaluation  of  Sl/j 

We  must  now  decide  how  to  evaluate  the  local  part  of  the  potential 

Siij,(x,t)  =  /  /  K(x,x' ,t  -  t')fi(x' ,t')dx' dt'  . 

Jt—s  Jr(t') 

Since  t  ~  t'  <  8,  the  kernel  K  is  sharply  peaked  at  x'  =  x.  This  rapid  decay  of  the  heat 
kernel  in  space  suggests  that  Scpt(x,t)  can  be  well-approximated  by  considering  only  the 
values  of  f.i  in  a  small  space-time  neighborhood  of  x  £  T(t)  as  8  — ►  0.  To  take  advantage 
of  this  locality,  we  expand  r(f')  and  n  in  Taylor  series,  and  construct  an  asymptotic 
approximation  to  Sl^l  in  powers  of  8. 

First,  assume  for  simplicity  that  T(t')  is  always  a  distance  >  d  from  the  boundary  of 
tiie  unit  box.  Then  K  can  be  approximated  by  the  free-space  kernel 

-|x-xT/4  (t-t1) 


as  t  —  t'  <  8  — +  0.  Indeed,  G  is  the  term  in 


,-||x-crX'-2k||2/.l(<-C) 


v  '  An(t  —  t1)  ^  ^  14 

1  >  keZJ  <r,=±l 

with  a,  =  1  and  k  =  0.  The  remaining  terms  are  decaying  exponentially 


\K(x,x',t  —  t')  -  G(x  —  x',t  -  t')l  =  0(— — ) 


for  x,x'  on  r(f').  Thus, 


ft  r  e-d*/s 

SLfi(x,t)  =  /  /  G(x  -  x',  t  -  ,t')dx' dt'  4-  0( — - — ) 

Jt-sJ  r(e)  8 
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as  8  — >  0.  This  last  expression  is  invariant  under  Euclidean  motions,  so  that  for  ease  of 
calculation  we  may  assume  that  x  =  0  and  that  the  tangent  line  to  T(t)  at  x  =  0  is  the 
j-axis.  We  also  assume  that  r(f')  extends  to  infinity,  incurring  an  approximation  error 
which  decays  like  0(e~1^6)  as  8  — >  0.  Then  T(t')  may  be  parametrized  by 

x  =  s,  y^y(sj') 

wtmre  y(0 ,t)  =  ys(Q,t)  =  0,  k  =  yss(0,f)  is  the  curvature  and  v  —  yt(0,t)  is  the  normal 
velocity  of  T(t)  at  s  =  0. 

After  these  notational  simplifications,  Siy(0,t)  assumes  the  form 

rt  fo o  e-s2 /4(t-t')e-y2(s,t')/4(t-t')  - - 

85  LJ.X - SITTej - +  , 

where  we  have  written  y(s,t')  in  place  of  y(x',t')  for  simplicity.  Note  that  s  is  equal  to 
the  arclength  at  s  =  0.  The  changes  of  variable  s  =  \J^(t -  t1)  ■  r  and  z  =  yf-i(t  —  t')  give 

1  fCO  ^  o  , - 

SLy( /  /  e~T\-y  (*T't~z  y(zr,t-  z7 /A)Jl  +  yVzrJ  -  z2/4)drdz. 

2t  Jo  J-oo  v 

The  sharp  peak  of  the  Gaussian  in  r  and  the  short  interval  of  integration  in  r  allow  us  to 
compute  S[^y  as  an  asymptotic  series  in  8.  After  Taylor  expansion  of  y  and  y  about  the 
point  (0,t),  we  obtain  to  lowest  order  in  6  the  result 

S^(0,i)  =  ^(0,()  +  O(«3/!). 

Retaining  one  more  order  in  6  and  doing  some  algebra  gives 

SLy{0,t)  =  ^(y(0,t)  +  ~(k  -v)2y(0,t)  -  |( yt  -  y„))  +  0(85/2) . 

Since  the  parameter  s  is  equal  to  arclength  at  s  =  0  and  the  curvature  and  normal 
velocity  are  invariant  under  Euclidean  motions,  the  preceding  formula  holds  more  generally. 
We  have  therefore  proven 

Lemma  3.3  Let  T(t)  and  y(x,t)  be  four  times  differentiable.  Then 

SLy(x(s, t),t)  =  ^  (l  +  -^(/c-w)2)/i(3,t)-|(Ait(5»0-A‘«(si0)  +0(8%),  (26) 

where  k  is  curvature,  v  is  normal  velocity  and  s  is  arclength  on  T(i). 

This  calculation  can  be  extended  to  higher  order  in  8  for  sufficiently  smooth  curves 
and  densities,  but  we  shall  carry  it  no  further.  The  approximation  of  Lemma  3.3  will 
suffice  for  the  purposes  of  this  paper:  it  can  clearly  be  evaluated  at  M  points  in  O(M) 
work,  completing  our  derivation  of  an  O(MN)  algorithm  for  evaluating  the  single  layer 
heat  potential  Sy. 
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3.3  Formal  Description  of  the  Algorithm 

In  this  section,  we  describe  the  fast  algorithm  in  a  more  procedural  form.  We  first  observe 
that  the  parameter  8  was  chosen  in  subsection  3.1  to  be  of  the  order  M-2^3+<,  independent 
of  the  time  step  At.  In  practice,  it  is  convenient  to  have  8  be  equal  to  an  integer  number 
of  time  steps.  We  therefore  set 


<5  = 


jV/-2/3+c 

At 


•  At. 


Algorithm 

Comment  [Choose  M,N,At  and  the  small  parameter  c  which  defines  8.  Set  l  =  -"|  and 

8  =  lAt.  Set  p  =  A/1/3.] 

Step  1. 

Comment  [Use  the  local  approximation  Spp  to  evaluate  the  single  layer  heat  potential  Sp  for  the  * 

first  l  time  steps.] 

do  n  =  1,2,...,/ 

do  m  =  1, 2, . . .,  M 

Evaluate  Sp(xm,nAt)  via  the  approximation  for  S^/i(xm,  nAt) 
given  by  equation  (26),  with  nAt  in  place  of  8  and  t. 

enddo 

enddo 


Step  2. 


Comment  [Initialize  Fourier  coefficients.] 

do  fci  =  1,  ...,p 
do  k2  =  l,...,p 

k:=  (fci,fc2) 

Ck(/At,£)  :=  0 

end  do 
end  do 


Step  3. 

Comment  [For  all  subsequci  *.  js,  update  the  Fourier  coefficients  and  calculate  the  history  part  ( 

Spp  at  the  points  ;.m  on  the  curve.  Then  add  the  local  approximation  Spp  at  each  xm  to  complete 
the  evaluation  of  the  single  'aye-  r  potential  Sp .] 
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do  n  —  l  1,1  +  2,  ■  ■  ■ ,  N 


(A)  do  =  1, p 

do  k2  =  1  ,-,P 
k  =  {k\,  k2) 

Evaluate  the  trigonometric  moments  Mk((n  -  l)Af,£)  and  Mk(nAt, &). 

Compute  the  update  Uk(nAt,  At,  6)  of  equation  (19) 
by  product  integration  (equations  22  -  24). 

Uk(nAt,  At,  6)  =  W0  •  Mk((n  -  l)Af,  6)  +  W\  •  Mk(nAt,6). 

Update  the  Fourier  coefficients  by  means  of  equation  (20): 

Ck(nAt,8)  =  e-*2ikl 2*‘Ck((n  -  l)At,6)  +  Uk(nAt,At,6) 
end  do 
end  do 

(B)  Evaluate  Fourier  series  at  p  equispaced  points  on  T(nAt). 

Extend  values  of  Fourier  series  to  all  M  points  xm  by  interpolation. 

(C)  do  m  =  1,2, . . . ,  M 

Evaluate  Scp(xm,nAt)  by  the  local  approximation  (26)  and  add  to  SF^(xm,n At). 

end  do 

end  do 

Remark  3.2.  Inspection  of  the  above  algorithm  shows  that  the  amount  of  work 
required  is  of  the  order  O(iVM),  assuming  fixed  degree  interpolation  in  arclength  is  used. 
If  Fourier  interpolation  is  used  in  Step  3  (B),  the  amount  of  work  required  is  of  the  order 
0(NM  logM). 

Remark  3.3.  As  noted  previously,  the  convergence  rate  of  the  scheme  is  dependent  on 
the  product  integration  scheme  used,  the  interpolation  scheme  used  and  the  order  of  the 
local  approximation.  The  scheme  described  in  the  text  is  second  order  in  /At  and  of  order  | 
in  8.  We  have  chosen  in  our  tests  to  use  Fourier  interpolation  in  arclength,  which  increases 
the  computational  complexity  by  a  factor  of  log  M,  but  effectively  removes  interpolation 
as  a  limiting  source  of  error. 

Remark  3.4.  In  the  algorithm  outlined  above,  we  end  up  computing  each  trigonomet¬ 
ric  moment  Mk(nAt,6 )  twice,  once  for  the  update  Uk(nAt,  At,  8)  and  once  for  Uk((n  + 
l)Af,  At,  8).  This  is  easily  avoided  by  an  appropriate  modification  to  the  program. 


4  Numerical  Examples 

The  algorithm  was  implemented  in  FORTRAN  and  tested  on  several  numerical  exam¬ 
ples.  We  evaluated  the  single  layer  heat  potential  of  a  cosine  density  p,{9)  =  cos(fc0)  of 
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wavenumber  k  on  a  stationary  circle 

T  :  (x  =  RcosO,y  =  R  sin  0,0  <  6  <  27r)  . 

We  computed  the  potential  at  time  T  =  1/2  with  a  sequence  of  numerical  parameters  in 
which  the  number  N  of  time  steps  and  the  number  M  of  points  on  the  circle  are  doubled 
at  each  stage,  p  and  the  number  Mp  of  points  used  for  integrating  over  the  curve  grow 
like  A/1/3,  and  6  decreases  like  A/-2/3.  Thus,  the  total  error  should  be  dominated  by  the 
65/2  =  M-s/3 

error  due  to  Taylor  expansion  in  the  local  approximation.  The  parameters 
used  and  corresponding  computational  times  are  shown  in  Table  1. 


Case 

N 

M 

6 

P 

B9i 

Time(Fast) 

Time(  Direct) 

1 

10 

mm 

.01 

10 

20 

9 

45 

2 

E9 

B9 

.0063 

13 

40 

12 

331 

3 

El 

80 

.004 

16 

50 

21 

4400 

4 

80 

160 

.0025 

20 

63 

56 

66024 

5 

160 

320 

.0016 

26 

80 

213 

(1.1  106) 

6 

320 

640 

.00099 

32 

102 

794 

(1.7  107) 

7 

640 

1280 

.00063 

40 

128 

3074 

(2.7  108) 

Table  1 

Table  of  parameters  for  Cases  1-7,  with  CPU  times  on  the  Multiflo  Trace  computer  at  Yale  University. 

In  cases  1-4,  the  direct  computation  times  were  estimated  by  evaluating  the  potential 
at  20  of  the  boundary  points  at  each  step,  using  the  trapezoidal  rule  in  space  and  second 
order  composite  product  integration  in  time.  In  cases  5-7,  direct  CPU  times  were  estimated 
by  extrapolation. 

The  accuracies  of  the  fast  and  direct  methods  are  compared  in  Table  2  for  wavenumber 
k  =  0.  The  error  reported  is  the  maximum  deviation  of  the  computed  potential  from  the 
exact  potential  over  20  points  on  the  curve.  Finally,  Table  3  presents  the  error  produced 
by  the  fast  algorithm  for  wavenumbers  k  =  0, 1,2  and  3. 

The  following  observations  can  be  made  from  Tables  1-3. 

1.  The  0(As5/3)  error  displayed  by  the  fast  algorithm  agrees  with  the  error  analysis 
of  section  3.2.  The  direct  algorithm  displays  0(Af2  +  As2)  error,  but  with  a  larger 
constant  of  proportionality,  which  accounts  for  its  poorer  accuracy  in  the  range  of 
parameters  tested.. 

2.  The  CPU  time  requirements  of  the  fast  algorithm  clearly  grow  only  like  NM. 

3.  By  the  time  N  =  640,  M  =  1280,  the  fast  algorithm  is  about  88,000  times  faster 
than  the  direct  method  would  have  been.  Case  7  required  about  50  min.  on  the 
Multiflo  Trace.  The  direct  calculation  would  have  taken  8-  years. 
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Case 

fast  algorithm 

direct  algorithm 

1 

WKEESBMM 

2 

3.1  10-5 

3 

9.2  10-6 

1.3  lcr4 

4 

2.7  10-6 

3.5  10-5 

5 

8.9  10-7 

8.5  1C)-6 

6 

2.9  10"7 

2.1  10-6 

7 

1.1  io-7 

5.1  10"7 

Table  2 

Comparison  of  absolute  errors  for  wavenumber  A:  =  0,  using  fast  and  direct  algorithms. 


Case 

o 

II 

k  =  1 

k  =  2 

■B9 

1 

1.2  10~4 

1.1  10"3 

2 

3.1  IO'5 

ppiig| 

3 

9.2  IO”6 

4 

2.7  IO"6 

4.5  10~6 

3.0  IO"5 

2.6  IO"4 

5 

8.9  1CT7 

1.4  IO'6 

9.8  10"s 

8.6  IO"5 

6 

2.9  10~7 

4.2  IO"7 

3.0  10-6 

2.6  IO"5 

7 

1.1  IO"7 

1.3  10~7 

9.5  IO'7 

8.5  IO'6 

Table  3 

Table  of  absolute  errors  produced  by  the  fast  algorithm  for  Cases  1-7,  with  wavenumbers  k  =  0,  1, 2,  3. 
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4.  Even  for  as  few  as  10  time  steps  and  20  boundary  points,  tbe  fast  algorithm  is  about 
five  times  faster  than  the  direct  calculation. 

5  Conclusions 

In  this  paper,  we  have  developed  a  fast  algorithm  for  evaluating  heat  potentials.  A  de¬ 
tailed  description  of  the  analysis  is  given  only  for  the  single  layer  potential  in  two  space 
dimensions,  but  the  outline  of  the  method  is  the  same  in  one  or  three  space  dimensions. 
The  extension  to  double  layer  potentials  is  straightforward. 

Our  algorithm  evaluates  a  heat  potential  at  MN  points,  using  density  values  at  MN 
points,  in  O(MN)  work.  The  direct  evaluation  requires  0(M2N2)  work,  so  the  fast  al¬ 
gorithm  achieves  a  dramatic  speedup:  with  1280  points  on  the  curve  and  640  time  steps, 
the  fast  algorithm  ran  88,000  times  faster  than  the  direct  calculation  would  have,  and 
produced  an  error  five  times  smaller. 
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